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I .  INTRODUCTIOH 

This  report  describes  the  improvement  and  testing  of  two  objective 
analysis  methods  developed  to  provide  initial  conditions  for  the  AFGL  Global 
Spectral  Model  (GSM).  A  previous  report  (Gerlach,  1983^)  described  the 
construction  of  the  two  methods.  Subsequent  to  that  report,  modifications 
were  made  to  improve  the  procedures,  and  the  modified  procedures  were  then 
used  in  conjunction  with  the  GSM  to  generate  a  series  of  forecasts  in  a  data 
assimilation  sequence.  The  first  procedure,  a  multivariate  optimum  inter¬ 
polation  (01)  procedure,  was  patterned  after  the  01  scheme  of  the  National 
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Meteorological  Center  (NMC)  as  described  by  Bergman  (1979)  ,  McPherson  et 
al.  (1979)^,  Kistler  and  Parrish  (1982)^,  and  Dey  and  Morone  (1983)^. 

The  Gerlach  report  described  the  design  of  the  procedure  in  some  detail,  so 
only  the  modifications  made  since  that  report  will  be  covered  in  this  paper. 
The  second  method,  an  attempt  to  fit  data  to  basis  functions,  was  detailed  by 
Halberstam  and  Tung  (1984)^.  Again,  only  subsequent  modifications  and 
testing  will  be  discussed  in  the  present  report. 

1.  Gerlach,  A.  M. ,  ed.,  1983:  Objective  Analysis  and  Prediction  Techniques 

-  1983.  AFGL-TR-83-0333 ,  Contract  F19628-82-C-0023 ,  Systems  and  Applied 

Sciences  Corporation,  ADA142441. 

2.  Bergman,  K.  H. ,  1979:  Multivariate  analysis  of  temperatures  and  winds 
using  optimum  interpolation.  Mon .  Wea .  Rev . .  107 .  1423-1444. 

3.  McPherson,  R.  D.,  K.  H.  Bergman,  R.  E.  Kistler,  G.  E.  Rasch,  and 

D.  S.  Gordon,  1979:  The  NMC  operational  global  data  assimilation  system. 

Mon .  Wea .  Rev . .  107 .  1445-1461. 

4.  Kistler,  R.  E.  and  D.  F.  Parrish,  1982:  Evolution  of  the  NMC  data 
assimilation  system:  September  1978-January  1982.  Mon .  Wea .  Rev . .  110. 
1335-1346. 

5.  Dey,  C.  H.  and  L.  L.  Morone,  1983:  Evolution  and  performance  of  the 
National  Meteorological  Center  Global  Data  Assimilation  System: 
January-December  1982.  Submitted  to  Mon .  Wea .  Rev . 

6.  Halberstam,  I.  M.  and  S.-L.  Tung,  1984:  Objective  analysis  using  Hough 
vectors  evaluated  at  irregularly  spaced  locations.  Mon .  Wea .  Rev . .  112. 
1804-1817. 


II.  OPTIMUM  INTERPOLATION  PROCEDURE 

A.  Recent  Modifications  of  AFGL  Statistical  Analysis  Programs  (ASAP) 

The  optimum  interpolation  procedure  developed  for  AFGL,  designated 
the  AFGL  Statistical  Analysis  Programs  (ASAP),  involves  a  terrain  surface 
pressure  analysis  followed  by  a  separate  upper  air  mass,  motion,  and 
moisture  field  analysis.  The  surface  pressure  analysis  uses  surface 
observations  of  pressure  and  winds  to  correct  the  forecast  of  surface 
pressure  by  the  GSM  on  the  model's  terrain  surface.  The  upper  air  code 
performs  corrections  of  height,  wind,  and  specific  humidity  forecasts  on 
model  o(=  layers  using  several  types  of  upper  air  observations. 

The  following  paragraphs  discuss  the  changes  made  in  the  two  codes  since 
the  Gerlach  (1983)  report. 

The  only  changes  in  the  surface  pressure  code  involve  the  so-called 
"buddy  check"  procedure.  Formerly,  this  procedure  was  performed  after  the 
observations  that  might  affect  a  particular  grid  point  were  located.  The 
observations  grouped  around  that  grid  point  were  checked  by  comparing  their 
observation-minus-forecast  residuals  with  a  tolerance  that  depends  on  the 
distance  between  the  observations  (see  Bergman's  Eq.  7.1).  This  comparison 
is  still  made,  and  the  flags  still  set  accordingly,  but  now  this  check  is 
performed  on  all  observations  before  the  optimum  interpolation  procedure 
begins.  The  observations  are  grouped  in  5®  latitude-longitude  boxes,  and 
within  each  box  the  buddy  check  is  performed.  If  an  observation  is 
excluded  in  this  process,  it  is  excluded  entirely  from  the  interpolation 
procedure  that  follows.  This  eliminates  the  possibility  of  an  observation 
being  used  in  the  analysis  of  one  grid  point  and  being  "flagged  out"  of  the 
analysis  of  a  neighboring  grid  point.  The  assignment  of  "toss"  flags  is  as 
it  was  in  the  Gerlach  (1983)  report,  but  an  assignment  of  a  "keep"  flag  to 
the  higher  quality  observation  (or  if  both  are  of  equal  quality,  to  both 
observations)  was  added  for  each  case  where  a  toss  flag  was  not  assigned. 
After  toss  and  keep  flags  are  assigned  for  all  observations  in  a  particular 
box,  the  toss  flags  are  removed  from  all  observations  with  two  or  more  keep 
flags  before  the  iterative  procedure  to  remove  the  observations  with  the 
most  toss  flags  begins.  In  this  way,  no  observation  with  two  or  more  keep 
flags  is  ever  eliminated.  This  follows  the  practice  of  NMC  as  given  by 
Kistler  and  Parrish  (1982). 
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The  question  of  the  necessity  of  a  vertical  interpolation  of  the 
first  guess  fields  to  the  new  a  pressures  as  defined  by  the  updated 
surface  pressure  remains  controversial.  Originally,  the  experimental 
design  for  this  project  called  for  testing  both  the  inclusion  and  exclusion 
of  such  a  vertical  interpolation  after  the  surface  pressure  analysis  and 
before  the  upper  air  analysis.  Unfortunately,  lack  of  time  and  a  shortage 
of  remaining  computer  funds  made  this  impractical,  so  the  vertical 
interpolation  was  retained  for  this  experiment. 

The  new  buddy  check  procedure  described  above  for  the  surface 
pressure  analysis  was  also  included  in  the  upper  air  analysis.  For  height, 
zonal  wind,  meridional  wind,  and  specific  humidity  (Z,  u,  v,  q),  the 
residuals  are  grouped  as  complete  observations  (that  is,  as  soundings)  on 
a  layers/levels  in  5“  latitude-longitude  boxes.  Then  box  by  box,  the  Z, 
u,  V,  and  q  residuals  are  considered  separately  by  layers/ levels ,  with  each 
of  the  four  variables  first  being  subjected  to  a  gross  error  check.  The 
error  limits  against  which  the  Z,  u,  v,  and  q  residuals  are  checked  are 
defined  as  four  times  the  forecast  error  standard  deviation  values  listed 
in  Table  1,  interpolated  linearly  in  In  p  to  the  corresponding  a 
layer/level  pressures.  Values  in  Table  1  for  Z,  u,  v  were  taken  from  Dey 
(1983)^,  and  values  for  q  up  through  300  mb  were  compiled  from  NMC 
forecast  statistics  for  March  1983  (Morone,  personal  communication). 

Values  for  q  above  300  mb  were  arrived  at  by  assuming  that  forecast  error 
decreases  with  decreasing  magnitude  of  q.  As  before,  Z  residuals  are 
located  at  the  a  layer  interfaces,  or  levels,  while  u,  v,  and  q  are 
located  in  the  o  layers.  Any  residual  exceeding  the  gross  error  limit  is 
eliminated.  Then  the  toss  and  keep  flags  are  assigned  to  the  remaining 
residuals  based  on  failing  or  passing  the  tolerance  check  of  pairs  of  like 
variables  at  the  layer/level  in  question.  Residuals  with  more  than  two 
keep  flags  are  relieved  of  their  toss  flags,  and  then  the  iterative 
procedure  to  remove  residuals  with  the  largest  number  of  toss  flags 
begins.  Each  succeeding  layer/level  is  checked  in  this  way  until  all  12 
levels  are  checked;  then  the  next  box  is  examined  in  the  same  way. 

7.  Dey,  C.  H. ,  1983:  The  NMC  optimum  interpolation  procedure.  Presented  at 
NMC  workshop  on  vector  processing  and  the  statistical  analysis  of 
meteorological  data.  Camp  Springs,  MD. 


Another  major  change  made  in  the  height-wind  multivariate  analysis  is 
the  use  of  the  same  residuals  to  calculate  a  correction  for  all  of  Z,  u,  and 
V  at  a  particular  grid  point  for  a  particular  a  layer.  The  proximity  of 
each  height  (on  levels)  and  wind  (on  layers)  residual  with  the  a  layer 
grid  point  to  be  corrected  is  calculated  using  the  product  of  the  non- 

dimensional  distance  functions  and  given  by 


"U  -  I-S, 

where  i  is  the  index  of  the  observation,  g  the  index  of  the  grid  point,  k.  = 

-12  -2  " 
1.96  X  10  m  ,k  =5.0, p  pressure  and  AS.  the  map  distance  between 

P  is 

the  observation  and  the  grid  point.  The  N  (a  variable,  currently  set  at  10) 

height  or  wind  observations  with  the  largest  value  of  are  then 

ig  ig 

chosen  for  the  calculation  of  the  correction  for  Z,  u,  v.  Thus,  since  both 
wind  components  are  taken  together,  as  many  as  2N  individual  residuals  may  be 
used  to  calculate  a  correction  for  each  of  z,  u,  and  v  at  the  a  layer  grid 
point.  Using  the  same  residuals  for  correcting  all  three  variables  has  the 
two-fold  advantage  of  enhancing  the  multivariate  nature  of  the  analysis  while 
saving  computational  time  since  only  one  correlation  matrix  is  formed  for  the 
correction  of  all  three  variables. 

Next,  the  actual  correlations  between  the  selected  residuals  and  the 

grid  point  a  layer  Z,  u,  v  are  calculated.  These  form  the  three  right- 

hand-  side  vectors  which,  when  combined  with  a  matrix  made  u^  of  correlations 

between  the  selected  residuals,  make  up  the  three  equation  sets  for  the 

corrections.  The  systems  of  equations  are  solved  for  the  interpolation 

weights  for  Z,  u,  v,  respectively,  using  the  Cholesky  method  for  solving  a 

0 

series  of  linear  equations  (Stobie,  1984  ).  Once  the  non-dimensional 
weights  forming  the  three  solution  sets  are  dimensionalized,  they  are  used 
along  with  the  corresponding  residuals  to  calculate  the  weighted  sum  which  is 
the  correction  to  be  applied  to  the  first  guess  values  of  Z,  u,  v, 

8.  Stobie,  J. ,  1984:  Cholesky  Method  for  Solving  a  Series  of  Linear 
Equations .  TSIN  Office  Note  84-1,  Air  Force  Global  Weather  Central,  Offutt 
AFB,  NE. 
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respectively,  at  that  a  layer  grid  point.  The  non-dimensionalized  weights 
are  used  with  the  respective  right-hand-side  vectors  of  residual-grid  point 
correlations  to  calculate  the  normalized  analysis  error  (see  Bergman's  Eq. 
2.13). 

Since  temperature  on  the  o  layers  is  the  mass  variable  in  the 
prognostic  model,  the  height  corrections  on  the  layers  have  to  be  converted 
to  temperature  corrections  on  the  layers.  First,  the  layer  height 
corrections  AZ  are  converted  to  mean  temperature  corrections  between  the 
13  layer  positions  of  height  (an  extra  height  correction  is  obtained  from 
the  01  height  analysis  for  a  layer  of  thickness  Ao^^  below  the  terrain 
surface)  using 

The  mean  temperature  corrections  are  assigned  to  the  ’’layer- layer”  pressures 
at  the  12  intermediate  points  between  the  13  a  layer  pressures.  The  a 
values  for  these  intermediate  layer-layer  positions  are  defined  such  that 
the  natural  logarithm  of  the  layer-layer  sigma  is  the  arithmetic  average  of 
the  natural  logarithm  of  the  sigma  values  of  the  surrounding  layers.  The 
layer- layer  mean  temperature  corrections  are  then  converted  to  a  layer 
temperature  corrections  using  the  Flattery  algorithm  (from  the  NMC  global 
spectral  model  preprocessing  code)  detailed  in  Appendix  A.  This  algorithm 
has  been  used  to  convert  interface  temperatures  to  layer  temperatures ,  and 
is  presented  in  this  way  in  Appendix  A.  The  application  of  this  technique 
to  the  problem  of  converting  layer-layer  temperature  corrections  to  layer 
temperature  corrections  is  carried  out  analogously,  and  in  our  case  K  =  13, 
because  of  the  use  of  the  subsurface  layer  to  avoid  extrapolation  in  the 
Flattery  procedure  in  the  lowest  model  layer.  By  calculating  a  correction 
for  the  subsurface  layer  in  the  upper  air  01  analysis,  the  residuals  are 
extrapolated  to  the  pressure  corresponding  to  using  the  vertical 
structure  functions  in  the  same  way  that  extrapolations  may  occur  at  a 
layers  above  the  highest  nearby  residuals.  A  similar  01  extrapolation  to  a 
layer  above  the  highest  a  layer  was  not  done  due  to  the  sparsity  of  data 
and  the  lesser  sensitivity  of  temperature  corrections  to  corresponding 
height  corrections  at  those  altitudes.  This  upper  level  extrapolation  by 
the  Flattery  routine  can  occasionally  produce  temperature  corrections  in  the 
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top  layer  that  may  not  be  completely  realistic,  but  by  extrapolating 
corrections  rather  than  full  temperature  values,  there  is  no  systematic 
warming  of  the  upper  layer  due  to  the  positive  lapse  rate  at  this  altitude. 
Corrections  may  have  positive  or  negative  lapse  rates,  so  in  some  places  the 
corrections  would  represent  large  temperature  increases  while  in  other 
places  large  temperature  decreases. 

McPherson  et  al.  (1979)  describes  the  procedure  of  allowing  the 
estimated  prediction  error,  used  to  normalize  the  observation  errors  (see 
Bergman’s  Eq.  2.10d),  to  evolve  over  succeeding  cycles  of  a  data 
assimilation.  The  estimated  prediction  error  for  one  cycle  is  obtained  from 
the  estimated  analysis  error  of  the  previous  cycle.  An  unfortunate  drawback 
of  such  a  scheme  is  that  over  time,  sharp  gradients  of  estimated  prediction 
error  develop  between  regions  of  high  data  density  and  regions  of  low  data 
density  (between  continents  and  oceans,  for  example).  This  violates  the 
assumption  of  horizontal  invariance  of  forecast  error  standard  deviation 

used  to  derive  structure  function  relationships  used  in  the  01  analysis. 

.  9 

An  alternative  approach  suggested  by  Lorenc  (1981)  is  to 

accumulate  statistics  of  estimated  analysis  error  over  long  periods  of 
forecast-analysis  cycles,  and  use  these  along  with  forecast  error  statistics 
to  generate  an  estimated  prediction  error  field  that  may  remain  constant  in 


time  in  ensuing  runs.  If  A  represents  an  analysis  value  of  variable  r 

at  grid  point  g  and  P  represents  a  corresponding  prediction  value,  and 

a  8^  ^  pa  1/2 

if  E  is  the  estimated  analysis  error  calculated  from  E  =  E*^  (c  ) 

gr  gr  gr  gr 

where  E^^  is  the  estimated  prediction  error  actually  used  at  the  grid 

point  and  is  the  non-dimensional  analysis  error  (defined  by  the 

right-hand-side  of  Bergman's  Eq.  2.13),  then  a  new  value  of  E^  can  be 

gr 

calculated  using 


(eP’)2 


(A  _  p  )2  +  (E®  )^ 
gr  gr  gr 


The  overbars  represent  averages  over  many  forecast-analysis  cycles,  enough 
to  include  seasonal  and  year-to-year  variations  in  the  forecasts. 


9.  Lorenc,  A.  C.,  1981:  A  global  three-dimensional  multivariate  statistical 
interpolation  scheme.  Mon.  Wea.  Rev. .  109.  701-721. 
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these  observation  o  layer  temperatures  are  used  to  calculate  observation 
a  level  heights  using  the  hydrostatic  equation.  This  check  was  removed 
for  the  calculations  of  the  GSM,  INIT,  and  ANAL  curves  on  Figs.  1  (a-c). 
Thus,  a  bad  height  observation  resulted  in  erroneous  a  layer  temperature, 
which  introduced  an  error  of  equal  magnitude  in  height  at  each  level  above 
when  the  hydrostatic  equation  was  used  to  get  observation  o  level 
heights.  The  check  of  derived  observation  a  layer  temperatures  was  put  in 
to  avoid  this  vertical  propagation  of  error  in  calculating  residuals  for  the 
analysis,  and  should  have  been  left  in  for  the  RMS  error  calculations.  In 
the  other  five  cases,  the  error  grew  with  height  without  any  sign  of  poor 
quality  data  in  the  sounding--apparently ,  a  legitimate  bad  fit  of 
observatic.i  with  field  value. 

In  contrast  to  RMS  errors  for  height,  similar  curves  for  vector  wind 

error  and  surface  pressure  (Figs.  1  (d-g))  do  appear  to  show  an  increase  of 

error  with  time.  For  the  winds,  the  rate  of  increase  of  error  seems  to 

increase  with  altitude,  as  is  seen  by  comparing  the  general  slopes  of  the 

three  ASAP-based  curves  in  Figs.  1  (d-f).  One  likely  explanation  for  this 

drift  away  from  the  observations  is  that  the  values  used  tor  estimated 

prediction  error  (in  this  case,  the  values  from  Table  1)  may  have  been  too 

small  for  a  12h  forecast.  Thus,  the  normalized  observation  error  c? 

ir 

would  have  been  too  large,  so  that  observations  would  have  had  insufficient 
influence  in  the  correction.  Since  the  values  in  Table  1  are  based  on 
comparing  6h  forecasts  with  observations,  they  should  have  been  augmented  to 
some  degree  to  serve  as  12h  forecast  error  values.  Since  a  12h  forecast  is 
two  sequential  6h  forecasts  with  no  correction  performed  in  between,  one  way 
to  augment  the  values  in  Table  1  would  be  to  use  forecast  error  growth  rates 
to  estimate  how  large  the  error  should  be  after  six  more  hours.  Such  values 
were  not  available  for  this  experiment,  although  estimates  could  be  derived 
from  Table  1  of  McPherson  et  al .  (1979).  Future  runs  to  experiment  with  the 
proper  specification  of  estimated  prediction  error  should  be  undertaken  to 
eliminate  the  drift  away  from  observations  if  a  12h  cycle  is  maintained. 
Changing  to  a  6h  cycle  may  very  well  eliminate  the  drift,  since  the  actual 
forecast  error  in  each  cycle  would  be  less.  In  the  long  run,  prediction 
error  statistics  could  be  generated  to  determine  the  values  that  should  be 
used  for  estimated  prediction  error  for  this  particular  forecast  model  and 
analysis  scheme. 


'' 


field  was  being  used  to  determine  which  observations  would  be  used  in  the 
comparisons.  For  these  RMS  errors,  only  conventional  upper  air  observa¬ 
tions  (Type  1  observations  in  the  FGGE  II-B  data  set,  consisting  of 
rawinsondes,  pilot  balloons,  dropsondes,  etc.)  were  used. 

In  all  cases,  the  analysis  error  (ANAL)  curve  lies  below  the 
forecast  error  (GSM)  curve,  showing  that  the  analysis  is  indeed  effective 
in  bringing  the  forecast  closer  to  the  observations,  as  would  be  expected. 
The  slow  component  error  (INIT)  curve  lies  somewhere  in  between  the  two  in 
all  cases  except  for  V(og) ,  where  for  the  most  part  it  lies  slightly 
below  the  ANAL  curve.  In  all  but  this  case,  the  initialization  removes  the 
fast  modes  from  the  analysis,  but  in  so  doing  partially  negates  the  effect 
of  the  analysis  in  bringing  the  forecast  field  closer  to  the  observations, 
as  is  expected.  Evidently,  in  the  case  of  the  upper  level  winds,  the 
analysis  has  a  smaller  effect  on  the  forecast  field  than  the  initialization 
has  in  restoring  mass-motion  balance. 

The  RMS  errors  for  neight  in  Figs.  1  (a-c)  for  the  three  ASAP  fields 
show  a  high  degree  of  irregularity  with  time,  but  do  not  seem  to  reveal  a 
trend  toward  increasing  error.  At  all  three  levels,  the  RMS  error  takes  a 
jump  at  1/11/79  OOZ,  but  the  experiment  was  not  carried  out  long  enough  to 
see  if  this  was  part  of  a  discernible  trend.  Given  the  lower  RMS  errors  on 
either  side  of  the  values  for  1/11/79  OOZ  and  12Z,  it  appears  that  this  was 
due  to  bad  fits  to  data  at  a  few  data  points  rather  than  a  larger  error 
overall.  This  was  confirmed  by  examining  the  observation-minus-GSM  values 
for  1/11/79  OOZ.  Height  differences  were  excessive  at  about  12  observation 
sites  for  Z(Og) ,  while  the  observation-rainus-F3A  differences  at  those 
same  locations  were  normal.  When  these  12  observations  were  removed  from 
the  RMS  calculation,  the  RMS  error  for  GSM  [7.(0^))  at  1/11/79  OOZ 
dropped  back  to  a  level  comparable  to  values  for  earlier  times.  Thus,  for 
height  at  least,  most  of  the  irregularity  in  the  curves  can  be  ascribed  to 
bad  fits  at  just  a  few  locations  where  the  F3A  fit  with  observations  was 
normal.  In  seven  of  the  12  cases,  the  bad  fit  was  traced  to  a  bad  height 
observation  value  at  levels  below  the  affected  level.  In  the  calculation 
of  observation-minus-field  values  for  the  01  analysis  and  the  F3A,  a  check 
is  made  of  the  observation  value  temperatures  as  they  are  calculated  for 
the  a  layers.  If  the  difference  between  an  "observation"  value  and  the 
field  value  of  temperature  is  greater  than  20°K,  it  is  excluded  before 
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P’ig.  1(f).  Same  as  Fig.  1(a)  except  V  at  Layer  8  (~225  mb) 
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ANALYSIS  TIMES  (GMT) 

Fig.  1(e).  Same  as  Fig.  1(a)  except  V  at  Layer  4  (~575  mb) 
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Fig.  1(c).  Same  as  Fig.  1(a)  except  Z  at  Level  9  (~200mb) 


7/00  8/00  9/00  10/00  11/00  12/00 
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Fig.  1(b).  Same  as  Fig.  1(a)  except  Z  at  Level  5  (~500mb) 


Fig.  1(a).  RMS  Error  between  Four  Indicated  Fields  and  Con 
ventional  Upper  Air  Observations  for  Z  at  Level  3  (~800mb) 
Using  F3A  as  Basis  for  Observation  Selection 


of  the  forecast  (G-GI)  with  the  fast  component  in  the  analysis  suggests 
that  only  a  small  portion  of  the  analysis  fast  component  can  be  attributed 
to  the  forecast  (background)  field.  Since  any  fast  component  information 
introduced  in  the  mass-motion  imbalances  in  the  analysis  is  removed  in  the 
initialization,  only  the  slow  component  portion  of  the  information  supplied 
to  the  forecast  field  by  the  analysis  is  of  interest.  Therefore,  it  is  of 
importance  to  know  the  magnitude  of  the  slow  component  introduced  in  the 
analysis  (I-GI),  where  I,  the  initialized  analysis,  is  the  slow  component 
of  the  analysis,  and  GI  is  the  slow  component  of  the  forecast  field  cor¬ 
rected  by  the  analysis.  For  each  variable,  a  comparison  of  (A-I)  -  (G-Gl) 
with  vI-GI)  shows  the  relative  magnitude  of  the  fast  component  contributed 
by  the  analysis  to  the  slow  component  contributed  by  the  analysis.  The 
smaller  this  ratio  (in  parentheses  for  average  of  first  eight  cases  at  end 
of  I-GI  row),  the  better  the  analysis  is  in  contributing  useful  information 
to  the  initial  conditions  for  the  next  forecast.  In  the  first  eight  cases 
of  this  experiment,  it  appears  that  on  the  average  the  analysis  fast  mode 
to  slow  mode  ratio  is  lower  for  heights  at  the  low  levels  and  lower  for 
winds  at  the  higher  levels.  This  would  seem  to  indicate  that  the  initial¬ 
ization  procedure  adjusts  the  winds  more  than  the  heights  to  achieve 
balance  at  low  elevations,  but  modifies  the  heights  more  than  the  winds  to 
obtain  mass-motion  balance  at  the  higher  levels.  By  comparison,  surface 
pressure  values  received  much  less  slow  mode  information  from  the  analysis 
in  proportion  to  fast  mode  information  than  did  upper  heights  and  winds. 

That  the  analysis  procedure  is  effective  in  pulling  the  forecasts 
closer  to  the  observations  is  evident  from  Figs.  1  (a-g) .  The  RMS  errors 
upon  which  these  graphs  are  based  were  calculated  by  comparing  the  FGGE 
III-A  (F3A) ,  the  forecast  (GSM),  initialized  analysis  (INIT),  and  the 
analysis  (ANAL)  fields  with  observations.  The  first  step  in  each  case  was 
to  calculate  the  observation-minus-field  differences  using  F3A,  then  invoke 
the  gross  error  check  and  buddy  check  in  the  same  way  that  they  are  used  in 
the  analyses.  The  observations  that  survived  these  two  checks  formed  the 
basis  for  the  calculations  of  observation-minus-field  differences  for  GSM, 
INIT,  and  ANAL.  Thus,  residuals  for  these  three  fields  were  calculated  at 
the  same  observation  locations  and  sigma  layers/ levels  for  which  the  F3A 
had  valid  values.  By  using  the  F3A  in  this  way,  it  was  believed  that  the 
comparisons  of  the  other  three  fields  would  be  more  valid  since  a  neutral 
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Table  2.  Root  Mean  Square  Differences  Between  (1)  GSM  Forecast  Field  and  Subsequent 
Analysis  (G-A) ,  (2)  Analyzed  Field  and  Initialized  Analysis  (A-I),  (3)  GSM 
Forecast  and  Initialized  Analysis  (G-I),  (4)  GSM  Forecast  and  Initialized 
Forecast  (G-GI),  and  (5)  Initialized  Analysis  and  Initialized  Forecast  (I-GI) 


0.24  0,22  0.24  0.24  0.24  0.25  0.24  0.25  —  0.24 

2.38  2.52  2.57  2,56  2.63  2.68  2.64  —  —  2.55  (0.58) 
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of  1/7/19  OOZ  was  chosen,  and  a  FGGE  III-A  analysis  for  that  date  and  time 
was  interpolated  from  mandatory  levels  to  12  sigma  levels  and  truncated  to 
the  rhomboidal  30  spectral  representation  following  the  procedure  given  by 
Gerlach  (1983).  These  spectral  values  of  relative  vorticity,  divergence, 
temperature,  and  specific  humioxLy  at  the  12  sigma  layers  as  well  as  terrain 
surface  pressure  were  subjected  to  the  two  iteration,  four  vertical  mode 
NMl.  Resulting  initialized  fields  were  then  used  as  initial  conditions  for 
a  12h  forecast  using  the  GSM.  This  forecast  acted  as  the  first  guess  in  the 
implementation  of  the  ASAP  01  analysis  codes,  using  FGGE  II-B  observations 
for  1/7/79  12Z.  The  analyzed  fields  were  then  initialized  using  the  NMI;  a 
12h  forecast  was  then  run,  followed  by  an  analysis  for  1/8/79  OOZ.  This 
pattern  was  repeated  for  a  total  of  ten  forecast-analysis-initialization 
cycles,  ending  with  an  initialization  of  the  analysis  valid  for  1/12/79 
OOZ.  The  following  paragraphs  describe  the  results  of  the  data  assimilation 
experiment.  The  codes  were  run  on  the  Air  Force  Weapons  Laboratory  Cray 
computer,  taking  about  35  minutes  of  computer  time  per  cycle,  of  which  about 
30  minutes  were  due  to  the  analysis. 

Table  2  is  a  display  of  the  root  mean  square  (RMS)  differences 
between  the  respective  pairs  of  the  three  fields  generated  in  each 
assimilation  cycle.  For  vector  winds,  the  RMS  difference  is  the  magnitude 
of  the  difference,  defined  by 


1 


(V. 


a. 
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1/2 


The  difference  G-A  represents  the  change  i.;jposed  on  the  forecast  field  by 
the  analysis,  A-I  is  the  fast  component  of  the  analysis  (due  primarily  to 
mass-motion  imbalances)  since  the  initialized  analysis  contains  only  the 
slow  mode  of  the  analysis,  and  G-I  is  the  overall  change  due  to  the  analysis 
and  initialized  field.  The  fourth  RMS  difference  shown,  G-GI,  is  the  fast 
component  of  the  forecast  field  since  the  Initialized  forecast  contains  only 
the  slow  component  of  the  forecast  field.  The  G-GI  RMS  values  are  shown 
only  for  the  first  eight  cases  because  these  calculations  were  performed 
later,  and  at  that  time  fields  for  the  last  two  cases  were  not  retrievable 
because  of  computer  hardware  problems.  A  comparison  of  the  fast  component 
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since  no  such  statistics  are  available  yet  for  the  AFGL  data 

assimilation  system,  values  of  forecast  error  standard  deviation  from  Table 

1  are  used  for  estimated  prediction  error  These  values  are  first 

linearly  interpolated  latitudinally  to  the  analysis  grid  latitudes  assuming 

that  the  values  given  in  the  table  correspond  to  latitudes  of  QO^S,  0®, 

20®N,  40®N,  and  90®N,  respectively,  for  each  of  five  columns  in  the  table 

for  Z,  u,  and  v.  Then  in  the  analysis  for  each  grid  point,  the  resulting 

values  are  interpolated  linearly  in  In  p  to  the  a  layer  grid  points  for 

each  layer.  These  values  are  used  as  an  estimate  of  f.  for  each 

gr  ir 

residual  i  of  type  r  used  in  the  correction  of  the  grid  point  g  in  Bergman's 
Eq .  2.10d.  That  is,  no  attempt  is  made  to  further  interpolate  estimated 
prediction  error  to  the  observation  site  to  provide  the  denominator  for  Eq. 
2.10d--the  grid  point  value  for  that  variable  is  used.  The  numerator  of  Eq. 
2.10d,  the  estimated  observation  error  E?^,  is  interpolated  vertically 
(linearly  in  iln  p  to  the  pressure  of  the  observation)  or  extracted 


directly  from  Dey  and  Morone's  (1983)  Table  3.  Thus  the  normalized 
observation  error  for  the  observation  i  of  type  r  is  given  by  e?  = 

E°  /E^  . 

ir  gr  ^ 

After  the  normalized  analysis  error  c  is  calculated  using 

Bergman's  Eq.  2.13,  the  dimensionalized  analysis  error  is  obtained  using 

3  P  3  1  /  2 

E  =  E  (c  )  for  each  variable  r  at  grid  point  g.  The  E  values 
gr  gr  gr  ore 

for  Z,  u,  V,  q,  along  with  the  corrections  calculated  in  the  analysis  and 


applied  to  the  forecast  value  P  to  form  the  analyzed  value  A  ,  are 

gr  '  gr’ 

the  two  terms  whose  squares  are  averaged  over  many  cycles  to  form  updated 

p  *  3 

values  for  E  .  For  this  reason,  the  E  values  and  corrections  are 
gr  gr 


stored  from  each  analysis. 


B.  Global  Data  Assimilation  Experiment  Using  ASAP  Analysis 


A  global  data  assimilation  experiment  was  conducted  using  the  ASAP 
analysis  codes  in  conjunction  with  the  GSM  and  the  normal  mode  initial¬ 
ization  (NMI).  Observations  and  the  starting  analysis  for  this  experiment 
were  extracted  from  the  First  Global  Atmospheric  Research  Program  (GARP) 
Global  Experiment  (FGGE)  data  tapes^^.  A  starting  date  and  time 


10.  Obtained  from  Department  of  the  Air  Force,  OL-A,  USAF  Environmental 
Technical  Applications  Center  (MAC),  Federal  Building,  Asheville,  NC  28801. 
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In  all  of  the  graphs  except  the  Pg^:^  plot,  where  the  F3A  field 
used  was  that  obtained  from  the  height  and  temperature  analyses  using  the 
GETPS  procedure  (see  Gerlach,  1983,  Appendix  B) ,  the  F3A  error  was 
substantially  less  than  the  error  associated  with  ANAL.  This  would  imply 
that  the  FGGE  analysis  achieved  a  better  fit  with  the  observations  than  did 
the  ASAP  analysis.  A  large  part  of  the  difference  in  error  level  between 
F3A  and  ANAL  could  be  explained  purely  on  the  basis  of  length  of  forecast 
in  the  cycles  between  corrections.  In  order  to  obtain  an  estimate  of  the 
forecast  error  difference  between  a  6h  and  12h  forecast  using  the  GSM,  a  6h 
forecast  was  run  starting  with  the  same  initialized  FGGE  III-A  field  at 
1/7/79  OOZ.  The  6h  forecast  was  compared  directly  with  observations 
without  the  buddy  and  gross  checks.  The  number  of  radiosondes  available  at 
that  time  was  about  one-third  the  number  available  at  the  12h  interval 
observation  times.  Values  of  error  are  plotted  on  Figs.  2  (a-g)  for  the  6h 
GSM  forecast.  In  all  cases  for  height,  the  difference  between  GSM(6)  and 
F3A  at  1/7/79  OOZ  is  less  than  half  of  the  difference  between  GSM  (12h)  at 
1/7/79  12Z  and  F3A  at  1/7/79  OOZ.  While  the  claim  could  not  be  made  that 
an  analysis  of  the  6h  forecast  would  pull  the  error  down  to  the  level  of 
the  F3A  error,  the  net  error  would  be  substantially  less  just  by  performing 
more  frequent  corrections. 

Another  explanation  for  the  difference  between  the  levels  of  RMS 
error  for  the  ANAL  and  F3A  is  that  by  using  F3A  as  a  basis  for  selecting 
the  observations  for  performing  this  error  analysis,  preference  is  being 
given  to  observations  that  fit  best  with  the  F3A  field.  A  different  subset 
of  the  observations  would  survive  the  gross  and  buddy  checks  if  the  ANAL 
were  used  as  a  basis.  To  test  this  possibility,  the  ANAL  field  was  used  as 
a  basis  for  selecting  observations  for  the  RMS  error  calculations,  and  F3A 
was  evaluated  using  that  set  of  observations.  The  results  are  plotted  in 
Figs.  2  (a-g).  For  Z,  the  F3A  and  ANAL  error  levels  reversed,  indicating 
an  approximately  equal  fit  with  observations  between  the  two  analyses.  The 
results  for  V  and  Pg^^  remained  virtually  the  same.  For  winds,  this 
would  imply  that  the  FGGE  assimilation  fit  the  observations  better  at  all 
times  regardless  of  the  subset  of  observations  used  as  a  basis  for  the  RMS 
error  calculations.  In  the  case  of  Pg^^.  the  FGGE  IIl-A  analysis  of 
surface  pressure  was  not  used  as  a  basis  of  comparison- -the  GETPS-deri ved 
values  were  used  since  this  is  the  field  used  in  the  GSM.  Naturally,  this 
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Fig.  2(a).  RMS  Error  between  Four  Indicated  Fields  and  Conven 
tional  Upper  Air  Observations  for  Z  at  Level  3  (~800  mb)  Using 
ANAL  as  Basis  for  Observation  Selection 


field  would  not  fit  as  well  with  observations  as  would  an  actual  analysis 
such  as  the  ANAL.  Finally,  the  large  difference  in  error  level  between 
ANAL  and  INIT  for  indicates  that  the  analyzed  field  did  not  maintain 

a  mass-motion  balance  well  in  correcting  the  pressure  forecast.  The 
correction  imposed  by  the  analysis  was  modified  by  the  initialization  in 
order  to  restore  balance  to  the  analyzed  surface  pressure  field. 

A  third  set  of  RMS  error  graphs  (Figs.  3  (a-g))  shows  the  results  of 
using  the  observations  common  to  both  F3A  and  ANAL  basis  sets  for  heights, 
winds,  and  specific  humidity.  Here,  the  fits  to  the  observations  for  F3A 
for  heights  are  only  slightly  better  than  those  for  ANAL,  whereas  the 
relationships  between  the  curves  for  winds  show  very  little  change.  While 
F3A  and  ANAL  cases  are  not  shown,  very  little  relative  change  occurred  in  q 
between  the  three  bases.  T'^is  series  of  graphs  demonstrates  that  the 
comparison  of  two  or  more  analysis  methods  as  to  how  they  fit  observations 
can  depend  very  much  on  what  set  of  observations  is  chosen  in  the 
comparison  and  on  what  variables  are  being  compared.  The  results  of  the 
comparison  of  the  fields  with  observations  averaged  over  the  ten  cases  are 
shown  in  Table  3.  Notice  that  the  number  of  observations  in  the  bases 
decreases  from  the  F3A  basis  to  the  combined  F3A,  ANAL  basis.  The  latter 
set  presumably  contains  only  those  observations  that  are  agreeable  to  both 
F3A  and  ANAL  fields;  thus  the  RMS  curves  in  Figs.  3  (a-f)  are  much  smoother 
than  for  the  other  two  bases. 

In  addition  to  a  comparison  of  a  sequence  of  analyses  from  ASAP  and 
FGGE  III-A,  48h  forecasts  based  on  the  respective  analyses  for  1/9/79  OOZ 
were  conducted.  The  resulting  forecast  fields,  designated  ASAP48  and 
F3A48,  respectively,  were  compared  to  1/11/79  OOZ  observations  using  the 
F3A  analysis,  ASAP  analysis,  and  combined  sets  for  1/11/79  OOZ  as  a  basis 
for  selecting  observations  to  be  used  in  the  RMS  error  calculations. 

Results  are  displayed  in  Table  4,  along  with  comparable  values  for  the 
initial  conditions  for  both  forecasts  (INIT0900,  initialized  ASAP  analysis, 
and  F3A  (init.),  initialized  FGGE  III-A  analysis].  The  purpose  of  this 
experiment  was  to  see  if  any  anomalies  were  created  by  the  ASAP  analysis 
that  would  cause  a  longer  forecast  to  show  a  radical  departure  from 
reality;  that  is,  to  see  if  a  longer  forecast  based  on  ASAP  would  be 
stable.  Although  it  would  have  been  desirable  to  conduct  several  such 
forecasts  for  several  different  dates  and  times,  fiscal  constraints  ruled 
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Fig.  3(a).  RMS  Error  between  Four  Indicated  Fields  and  Conven 
tional  Upper  Air  Observations  for  Z  at  Level  3  (~800  mb)  Using 
Observations  Common  to  Both  F3A  and  ANAL  Basis  Observations 


C.  The  Iterative  and  Non-Iterative  Cycles 


Once  the  coefficients  ^  are  determined,  it  is  possible  to  aprly 
them  to  evaluating  the  variables  at  grid  locations.  At  this  point,  the 
procedure  differs  little  from  other  analyses.  To  perform  an  initialization, 
the  variables  must  be  expressed  in  terms  of  divergence,  vorticity,  tempera¬ 
ture,  surface  pressure,  and  moisture.  Velocities  expanded  at  grid  points  can 
be  rewritten  in  terms  of  divergence  and  vorticity.  The  compound  function  P 
can  be  separated  into  a  surface  pressure  and  height  component,  while  moisture 
is  provided  independently  from  the  FGGE  III-A  analyses.  Heights  can  be 
transformed  into  temperatures  using  a  hydrostatic  relationship,  but  the  pro¬ 
cedure  is  by  no  means  straightforward.  In  the  NMC  initialization  procedure, 
Sela’s  (1980)^^  matrix  equation  relating  temperatures  to  heights  is  in¬ 
verted  to  solve  for  temperatures  from  heights.  Yang  (1982)^^  has  pointed 
out,  however,  that  this  procedure  is  faulty  because  the  matrix  inversion 
results  in  unrealistic  temperature  profiles.  NMC  is  not  affected  much  by 
this  shortcoming  because  during  the  course  of  initialization,  the  compound 
function  P  changes  very  little,  and  it  is  only  the  change  that  is  computed  by 
the  faulty  matrix  inversion.  But  in  this  study,  the  actual  height  values  are 
transformed  back  to  temperatures  before  initialization,  and  this  presents 
problems.  As  in  Section  II  of  this  report,  the  Flattery  least-squares  fit, 
as  detailed  in  Appendix  A,  is  invoked  to  determine  temperatures  from  heights 
at  the  layers.  This  method  had  its  drawbacks,  as  well,  mentioned  in  the 
earlier  section  of  this  report.  These  faults  can  become  quite  serious  as 
will  be  explained  later. 

Once  all  the  required  variables  are  defined  properly  at  all  layers, 
the  initialization  procedure  can  be  performed.  Although  the  analysis  employs 
eight  vertical  modes,  the  initialization  limits  the  modes  to  four.  This  is 
to  prevent  divergence  of  the  Machenhauer  non-linear  iteration  scheme.  Once 
the  initialization  routine  has  balanced  the  mass  and  motion  field,  there  are 
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17.  Yang,  C.H.;  1982;  On  the  solution  of  the  hydrostatic  relation  in  the 
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coefficient  in  the  expansion.  The  second  coefficient  is  found  by  removing 
the  contribution  of  the  first  part  and  repeating  the  process;  i.e.. 


If  the  G  were  truly  orthogonal  to  one  another,  Eq .  (3)  would  reduce  to 
n ,  r 

the  simple  equation  governing  the  evaluation  of  coefficients  of  orthogonal 
expansions.  The  departure  from  orthogonality  will  determine  both  the 
accuracy  and  validity  of  Eq .  (3)  as  a  means  of  determining  the  coefficients 
b 

n.r 

As  mentioned  by  HT,  the  problem  of  ordering  of  the  vectors  ©  j. 
not  trivial.  It  is  possible,  as  we  attempted  to  do  in  this  study,  to  order 
the  vectors  by  their  contribution  to  the  data  vector,  i.e.,  by  the  magnitude 
of  This  procedure  effectively  doubled  the  analysis  time  on 

AKWL's  Cray  computer,  and  the  resulting  improvement  was  minor.  Because  we 
are  dealing  with  residuals,  rather  than  absolute  data,  it  was  not  possible 
to  predict  which  vector  would  have  the  greatest  impact.  In  fact,  the  vector 
O®,  which  is  a  constant  and  would  result  in  the  mean  of  each  variable, 
would  certainly  have  been  the  greatest  contributor  to  a  field  of  variables 
such  as  heights  or  temperatures.  In  dealing  with  residuals,  however,  there 
is  no  guarantee  that  the  mean  is  necessarily  greater  than  any  of  the  other 
amplitudes  of  the  expansion.  It  may  even  be  possible  that  vectors  from 
internal  modes  with  smaller  effective  depths  would  result  in  larger  ampli¬ 
tudes  than  the  external  mode  or  some  of  the  larger  internal  modes.  But  to 
order  all  I(L+1)(M+1)  vectors  would  prove  to  be  computationally  impracti¬ 
cal.  It  is  thus  necessary  to  separate  vertical  modes  and  to  deal  with  only 
the  (L+1)(M+1)  vectors  corresponding  to  each  vertical  mode.  It  was  found 

that  in  general  F  •  0.  >>  0  •  ©  for  any  value  of  s  and  k  as  long 

It  KjI/  S«P  iCjP 

as  s4k.  But  this  inequality  depends  greatly  upon  the  geographical  distri¬ 
bution  of  data  and  cannot  be  generalized  for  all  cases. 


In  theory,  I  could  be  as  large  as  the  number  of  a  levels  in  the  model,  in 
this  case,  12.  In  practice,  however,  only  the  first  eight  vertical  modes 
were  retained  in  the  expansion.  In  case  of  observations  with  missing 
levels,  the  residuals  at  those  levels  were  set  to  zero. 

After  the  horizontally  dependent  coefficients  were  determined  for 

each  vertical  mode,  they  were  expanded  in  horizontal  nmf  by  the  method 

described  by  HT.  This  meant  evaluating  the  nmf  at  each  observation  point 

where  the  residuals  were  defined.  As  with  HT,  this  included  the  corners  of 

grid  boxes  that  lacked  observations,  which  were  assigned  residuals  of  zero. 

Unlike  HT,  however,  this  study  included  single  observations  of  either 

heights  or  velocities,  which  were  fit  to  either  U,  V,  or  P  in  the  hope  that 

the  final  analysis  would  create  a  balance  between  the  heights  and  velocities 

even  in  regions  where  one  or  the  other  was  missing.  The  projected  residuals 

are  then  fit  sequentially  to  the  nmf  for  each  mode,  and  each  successive 

contribution  subtracted  from  the  original  projected  data  vector.  In  other 

words,  if  F  represents  a  vector  containing  the  coefficients  {a(d>. ,  X.)} 
r  1  j 

for  a  particular  vertical  mode  r  at  all  latitudes  and  longitudes  where 
residuals  of  height  and  velocity  are  found,  and  we  are  also  given  a  set  of 
vectors 


K.r)  ■ 


for  m  =  0,...,M  zonal  wavenumbers  and  9.  =  0,...,L  frequencies  evaluated  at 
the  respective  latitudes  and  longitudes  where  the  corresponding  u,  v,  and  P 
residuals  occur,  then  we  can  fit  the  vector  to  the  set  {0^  in  a 
sequential  manner.  We  first  order  the  set  {  0^  }  in  some  logical 
fashion,  allowing  us  to  substitute  one  index,  say  n,  for  the  H  and  m.  The 
index  n  will  then  range  from  1  to  (L+1)(M+1).  The  first  coefficient  is 
determined  by  ignoring  the  nonorthogonality  of  the  vectors  in  the  set  defined 
over  the  observation  sites  and  allowing 


b,  =  ( F  •  0,  )(6,  •  6,  t  where  b,  is  now  the  first 

l,r  r  l,r  l,r  l,r  l,r 
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where  f ,  d"*  .  ,  p"*  ,  are  normalized  spectral  coefficients  of  vorticity, 
n,a  n.l  -n.J. 

divergence,  and  the  compound  variable,  respectively,  for  zonal  wavenumber  m, 

meridional  wavenumber  n,  and  frequency  index  I,  is  the  Legendre  function, 

U^,  V^,  and  P^  are  the  normal  mode  functions  (nmf)  for  the  zonal  component  of 

velocity,  the  meridional  component,  and  the  compound  variable,  respectively. 

o'"  =  (n^-  (4n^-  1)  <J>  is  latitude,  and  i  E  /-I.  At  the  poles, 

n 

and  will  be  zero,  except  for  m  =  1,  when  P^(cos  <J))  ^  approaches 
*■  *■  1/2  — 

±  [n(n  +  1)  (2n  +  1)]  /2v'2  as  4>  approaches  ±  ir/2.  For  m  =  0,  the 

Rossby  waves  have  to  be  derived  separately  under  the  constraints  that 

V?  =  0  and  the  normalization  requirement  that  Z]  |C^  +  |P^  .1^  =  1.  A 

a  ^  n  n,a  -n,a 

set  of  solutions  is  then  obtained  by  means  of  Kasahara's  (1978)^^  procedure 
of  arbitrarily  allowing  one  coefficient  to  equal  1  and  computing  the  others. 
This  leads  to  a  sequence  of  vectors  which  can  then  be  orthonormal ized  with  a 
Gramm-Schmidt  procedure,  and  the  resulting  coefficients  are  combined  to  pro¬ 
duce  and  P^for  the  Rossby  modes. 

In  order  to  avoid  considerable  computational  efforts  in  finding  the 
values  of  the  nmf  for  every  observation  point,  we  calculated  and  stored  their 
value  for  each  degree  of  latitude  and  interpolated  whenever  necessary.  For 
high  wavenumbers  this  may  engender  some  error,  especially  near  the  equator, 
but  hardly  ever  greater  than  1  percent,  which  is  acceptable  for  this  study. 

A  suggestion  will  be  offered  later  as  to  how  to  reduce  this  error. 

To  simplify  the  three-dimensional  aspect  of  the  analysis  problem,  all 
observations  were  interpolated  to  model  o  layers  by  means  of  an  01  method 
described  in  Section  II  (pp.  9-10).  The  first  guess  field  was  then  fonried  by 
spectral  expansion  of  the  forecast  at  the  observation  sites  and  the  residuals 
computed.  The  resulting  field  of  residuals  was  then  projected  on  the  model 
vertical  modes  {Z^^Co)},  i  =  l . I.  Because  the  vertical  modes  are  orthonor¬ 

mal  with  respect  to  the  model  o  structure,  the  computation  of  the  projec¬ 
tions  was  fairly  straightforward  so  that  any  variable  A(4>,X,o),  where  X 

is  longitude,  defined  at  a  particular  location  <J).  ,  X.,  o,  ,  was  represented  as 

1  J  K 
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optimum  interpolation  (01)  but  was  hindered  by  the  computational  requirements 
for  his  recommendations. 

By  employing  the  normal  mode  functions  directly  in  the  analysis  pro¬ 
cess,  we  hope  to  merge  some  of  the  features  of  initialization  with  objective 

analysis.  In  this  study,  the  normal  mode  functions  were  the  same  as  those 

13 

generated  by  Ballish  (1980)  for  NMC's  spectral  model.  These  functions  are 

determined  by  linearizing  the  tendency  equations  for  vorticity,  divergence, 

heights,  and  surface  pressure.  The  resulting  solutions  supply  eigenvalues 

corresponding  to  effective  heights  in  the  atmosphere  and  eigenvectors  of 

spectral  coefficients  for  vorticity,  divergence,  and  a  compound  variable 

combining  surface  pressure  and  geopotential.  The  eigenvectors  are  separated 

by  zonal  wavenumber  and  frequency.  The  frequencies  can  be  classified  as 

Rossby,  or  rotational,  waves,  eastward-propagating  gravity  waves,  or 

westward-propagating  waves.  For  this  study,  as  well  as  for  NMC's  purposes, 

gravity  waves  with  a  period  of  48  h  or  more  are  kept  along  with  all  Rossby 

frequencies.  All  others  are  discarded.  As  with  the  Hough  functions  de- 

14 

scribed  by  Kasahara  (1976)  ,  the  coefficients  can  be  used  to  create  merid¬ 

ional  normal  mode  functions  for  velocity  and  the  compound  variable.  The 
expansions  are  of  the  following  form: 


a. 
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global  model  rather  than  Hough  functions,  which  are  normal  modes  of  the 
Laplace  tidal  equations.  The  data  were  projected  onto  eight  of  a  possible 
12  vertical  modes,  then  projected  along  the  horizontal  normal  modes.  There 
were  two  sets  of  experiments.  In  the  first  set,  an  initialization  was 
performed  after  the  analysis  and  forecasts  generated  from  the  initialized 
field.  In  the  second  set,  the  initialized  field  replaced  the  forecast  field 
as  the  first  guess,  and  a  second  analysis  was  produced  based  on  the  observed 
data  and  the  initialized  field.  This  procedure  was  repeated  once  more  be- 
fore  a  final  initialized  field  was  readied  for  the  forecast  model.  This 
iterative  procedure  was  designed  to  create  a  field  which  would  reflect  both 
the  data  and  the  large  scale  motion  of  the  atmosphere.  In  mathematical 
terminology,  the  iteration  was  meant  to  converge  on  the  intersection  between 
the  "slow”  and  "data"  manifolds. 

B.  The  Fitting  Functions 

During  the  course  of  preparation  for  a  numerical  weather  forecast, 
data  are  moved  from  their  observation  locations  usually  to  a  fixed  grid  and 
are  balanced  for  the  forecast  model  by  initialization.  The  interpolation  of 
data  to  a  fixed  grid  is  a  process  referred  to  by  meteorologists  as  an  "an¬ 
alysis."  The  balancing  of  the  resulting  fields  is  termed  "initialization," 
indicating  that  the  field  of  variables  is  adjusted  to  serve  as  initial  con¬ 
ditions  for  the  model.  In  recent  years,  the  initialization  process  has 
centered  about  model  normal  modes,  which  are  the  solutions  to  the  linearized 
model  equations.  To  balance  the  tendencies  of  the  variables,  the  non-linear 
terms  are  also  included  in  some  fashion.  The  entire  process  is  labeled 
"non-linear  normal  mode  initialization"  (NLNMI).  It  is  very  possible  that 
the  analysis  procedure  may  work  at  cross  odds  to  the  initialization  require¬ 
ments  unless  the  two  are  somehow  combined.  Interpolation  to  specific  grid 
points  alters  the  nature  of  the  data  being  interpolated.  When  initializa¬ 
tion  is  performed,  the  data  are  again  modified,  raising  doubts  as  to  the 
fidelity  of  the  final  product  to  the  observations.  It  would  be  advantageous 
to  combine  the  processes  of  analysis  and  initialization  so  that  the  final 

product  does  not  depart  radically  from  the  original  observations.  In  fact, 
12 

Phillips  (1982)  believed  such  an  integration  was  possible  in  the  case  of 

12.  Phillips,  N.A. ,  1982:  On  the  completeness  of  multi  variate  optimum 
interpolation  for  large-scale  meteorological  analysis.  Mon .  Wea .  Rev . .  110, 
1319-1334. 
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111.  ANALYSIS  WITH  DISCRETE  NORMAL  MODE  FUNCTIONS 


A.  Introduction 

A  prescription  for  combining  objective  analysis  with  normal  mode  ini¬ 
tialization  was  offered  by  Halberstam  and  Tung  (1984)^  (HT,  hereafter)  who 
demonstrated  that  Hough  functions  evaluated  at  observation  sites  can  be  fit 
to  height  and  wind  observations  to  produce  an  analysis  that  can  be  con¬ 
sidered  "pre-initialized . "  In  order  to  make  the  fitting  procedure  economi- 
ically  feasible,  it  was  necessary  to  compromise  on  a  least-squares  fit.  It 
was  shown  that  when  a  procedure  such  as  the  one  outlined  by  Holmstrom 

(1963)^^  was  substituted  for  a  least-squares  fit,  the  resulting  error  is 

•t 

smaller  than  a  least-squares  fit  in  areas  of  few  observations.  Holmstrom* s 
procedure  is  a  sequential  computation  of  the  coefficients  of  the  expansion 
by  subtracting  from  the  original  observation  vector  the  fit  produced  by  each 
vector  from  the  set  of  basis  function  vectors.  As  can  be  anticipated,  the 
geographical  distribution  of  data  plays  an  important  role  in  the  size  of  the 
analysis  error.  This  makes  precise  estimation  of  the  error  a  difficult 
task.  It  was  shown  by  HT  that  the  error  is  related  to  the  non-orthogonality 
of  the  vectors;  namely,  the  degree  to  which  the  dot  products  of  the  vectors 
depart  from  zero.  Several  comparisons  were  made  by  HT  using  a  set  of  Hough 
functions  from  the  external  mode  only,  with  500  mb  radiosonde  data  from  FGGE 
II-B.  Their  results  showed  that  although  rms  differences  between  analysis 
and  data  were  quite  large,  the  sequential  method  did  control  the  errors.  An 
iterative  procedure,  where  the  resulting  analysis  was  substituted  as  the 
first  guess  field  in  the  calculation  of  residuals,  resulted  in  smaller 
errors  between  the  observations  and  the  final  analysis. 

In  this  study  several  modifications  were  made;  each  will  be  explained 
in  detail.  The  modifications  were  geared  to  a  global  analysis  based  on  all 
FGGE  II-B  data  which  are  not  rejected  by  a  checking  procedure.  This  would 
include  satellite-derived  heights  or  single  level  aircraft  data.  The  basis 
functions  were  chosen  to  be  the  normal  mode  functions  of  the  AFGL  spectral, 

»• 

11.  Holmstrom,  I.,  1963:  On  a  method  for  parametric  representation  of  the 
state  of  the  atmosphere.  Tellus.  15 .  127-149. 
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C.  Conclusions 


Though  hampered  by  the  objectivity  problem  of  using  observations 
whose  error  levels  are  not  accurately  known  as  a  basis  for  verification,  it 
appears  that  the  01  analysis  procedure  described  in  this  and  two  preceding 
reports  works  reasonably  well  in  producing  global  forecasts  in  conjunction 
with  the  NMI  and  GSM.  Although  the  10  cycle  data  assimilation  procedure 
showed  some  increase  of  error  with  time  in  the  sequential  forecasts  for 
winds  and  pressure,  this  may  be  due  to  an  underestimate  of  estimated 
prediction  error  for  the  12h  forecasts  used  in  the  assimilation 
experiment.  For  further  study,  a  6h  forecast  length  is  recommended,  even 
though  fewer  conventional  observations  are  available  at  06Z  and  18Z.  In 
addition,  it  is  recommended  that  observations  with  high  quality  control 
index  values  (poor  quality)  from  the  F(X;e  II-B  data  set  not  be  used  with 
this  procedure,  since  their  inclusion  (especially  in  conventional  upper  air 
observations)  often  results  in  erroneous  interpolations  to  a  layers/ 
levels  when  use  of  a  correct  observation  just  above  or  below  the  erroneous 
one  would  have  led  to  an  acceptable  interpolated  value.  Finally,  while 
results  for  moisture  were  not  discussed  in  this  report  because  of  the 
simple,  univariate  nature  of  the  moisture  analysis,  the  results  of  the 
moisture  analyses  in  this  study  will  be  used  as  a  standard  against  which  an 
improved  01  based  global  moisture  analysis  method  will  be  compared  in  a 
follow-on  study. 
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methodologies.  Fortunately,  Dey  and  Morone  (1983)  show  curves  for  RMS 
error  against  observations  for  both  initial  conditions  and  6h  forecasts  for 
1  1/2  years  of  assimilation  performances  of  the  NMC  Global  Data  Assimila¬ 
tion  System  (GDAS).  These  curves  were  qualitatively  averaged  over  their 
duration  to  obtain  the  values  for  the  GDAS  shown  in  Table  5.  The 


Table  5.  Qualitative  Comparison  of  ASAP^l^  with  GDAS^^)  rmse  Against 
Observations . 


Z  (SOOmb) 

GDAS  (SOOmb) 
ASAP  (-SOOmb) 

V  (2S0mb) 

GDAS  (2S0mb) 
ASAP  (~22Smb) 


GSM  IMIT  DIFF 

24  (m)  18  6 

35.2  23.7  11. S 

8.6  (m  s-1)  6.4  2.2 

11.3  9.7  1.6 


(1)  Values  based  on  average  of  663  radiosondes  over  the  globe  for  ten 

cycles:  1/7/79  to  1/12/79  (12h  forecasts) 

(2)  Values  based  on  102  northern  hemisphere  radiosondes  for  6h  forecast 

cycles:  7/81  -  12/82  (Dey  and  Morone,  1983) 


corresponding  ASAP  values  (based  on  the  ten  12h  cycles)  averaged  over  the 
time  period  and  based  on  the  intersection  set  of  observations  (F3A  +  ANAL 
basis)  were  taken  from  Table  3.  The  absolute  values  for  both  the  forecast 
(GSM,  6h  for  GDAS,  12h  for  ASAP)  and  the  corresponding  initialized  analysis 
resulting  from  that  forecast  field  (IMIT)  are  larger  for  the  ASAP  system  in 
all  cases.  But  considering  only  the  net  change  imposed  by  the  analysis- 
initialization  processes  on  the  forecast  (the  DIFF  column)  shows  that,  at 
least  qualitatively,  the  effectiveness  of  the  ASAP  system  is  competitive 
with  that  of  the  GDAS. 


out  that  possibility.  Thus,  while  it  is  somewhat  risky  to  draw  sweeping 
conclusions  from  one  comparison,  it  appears  that  the  ASAP  based  forecast 
has  RMS  errors  that  are  of  the  same  order  of  magnitude  as  F3A  based 
forecast  BIMS  errors. 

In  all  cases,  the  F3A48  forecast  verifies  better  against 
observations  than  does  the  ASAP48  forecast.  However,  since  in  most  cases 
the  initial  conditions  also  verify  more  closely  in  the  F3A  case,  a  more 
accurate  measure  of  forecast  performance  is  the  error  growth  during  the 
course  of  the  48h  forecast  in  each  case.  By  comparing  the  difference 
between  F3A48  and  F3A  (init)  with  the  difference  between  ASAP48  and  ASAP 
(init),  we  see  that  the  error  growth  for  heights  was  less  for  ASAP  only  in 
the  Z  (o^)  "F3A  only"  case,  but  in  all  cases  for  the  highest  level 
winds.  The  two  forecasts  are  about  even  in  error  growth  in  predicting 
surface  pressure.  Thus,  in  this  case,  the  ASAP  forecast  does  not  appear  to 
be  the  victim  of  any  anomalous  error  that  would  have  resulted  in  a  grossly 
large  verification  error  at  48  hours. 

As  in  the  case  with  the  12h  assimilation  cycles,  the  fit  of  fields 
to  observations  very  much  depends  on  the  set  of  observations  chosen.  This 
48h  experiment  shows  that  the  forecast  verifies  best  in  almost  all  cases 
against  the  most  restrictive  set  of  observations,  the  intersection  of  the 
F3A  and  ANAL  based  sets,  which  in  all  cases  contains  the  fewest  observa¬ 
tions.  Once  again,  the  difference  between  statistics  for  three  different 
sets  of  basis  observations,  especially  for  Z,  indicates  that  results  can 
depend  on  which  observations  are  selected  for  the  verification  of  forecasts 
and  analyses.  This  results  from  the  fact  that  the  errors  associated  with 
observations  reduce  the  objectivity  of  the  verifications,  since  no  accurate 
objective  measure  of  observation  error  is  known.  No  two  methods  of 
throwing  out  bad  observations,  whether  manually  or  automatically,  will 
result  in  the  same  set  of  observations  against  which  the  relative 
"goodness"  of  two  or  more  forecasts/analyses  can  be  compared.  All  that  can 
be  said  in  this  case  is  that  the  RMS  errors  associated  with  the  two 
forecasts  are  of  the  same  order  of  magnitude. 

Finally,  as  a  basis  of  reference  for  a  newly  developed  assimilation 
system,  it  is  of  interest  to  establish  at  least  a  qualitative  comparison 
with  a  well-established,  operational  assimilation  system  to  see  if  the 
level  of  performance  of  the  system  is  competitive  with  established 
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Table  4.  RMS  Errors  for  48h  Forecast  Experiment 


Variable 

Experiment 

F3A(init) 

ASAP ( ini t) 

F3A48 

ASAP48 

No .  Obs 

Z(m) 

F3A  only 

17.81 

19.94 

28.22 

33.10 

654 

^3 

ASAP  only 

17.61 

16.79 

26.45 

31.48 

642 

F3A  +  ASAP 

12.83 

15.24 

23.89 

28.57 

639 

F3A  only 

25.96 

37.75 

70.01 

78.02 

642 

ASAP  only 

24.85 

25.67 

68.65 

76.29 

636 

F3A  +  ASAP 

21.55 

24.26 

58.69 

66.92 

625 

F3A  only 

43.28 

62.43 

105.73 

126.17 

614 

^9 

ASAP  only 

41.90 

54.13 

103.31 

123.85 

616 

F3A  +  ASAP 

40.26 

53.59 

95.95 

116.70 

603 

V(m  s“^) 

F3A  only 

4.38 

4.75 

8.08 

9,13 

675 

ASAP  only 

4.46 

4.51 

7.98 

8.98 

665 

F3A  +  ASAP 

4.35 

4.49 

7.94 

8.97 

652 

F3A  only 

4.62 

6.37 

8.92 

9.67 

631 

‘’a 

ASAP  only 

4.52 

5.98 

8.94 

9.71 

603 

F3A  +  ASAP 

4.45 

5.94 

8.79 

9.52 

593 

F3A  only 

7.04 

10.67 

14.73 

14.24 

582 

‘*8 

ASAP  only 

6.90 

9.32 

14.31 

13.82 

565 

F3A  +  ASAP 

6.74 

9.27 

14.29 

13,83 

561 

F3A  only 

3.34 

3.30 

4.95 

6.09 

2736 

^sf  c 

ASAP  only 

3.56 

3.17 

4.87 

6.05 

2715 

F3A  +  ASAP 

— 

— 

— 

— 

— 
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Table  3.  RMS  Differences  between  Fields  and  U.A.  Observations 


(Averaged  over  10  12-hour  Assimilation  Cycles) 


Obs.  Allowed 

Obs . 

Allowed 

Obs.  Allowed 

by  F3A 

by  ANAL 

by  F3A 

+  ANAL 

RMS  No .  Obs . 

RMS 

No.  Obs. 

RMS 

No.  Obs 

Z3=3  C~ 

800  mb) 

GSM 

22.67(m)  690 

20.41 

683 

19.38 

679 

INIT 

19.63 

16.84 

15.37 

ANAL 

18.83 

15.27 

15.16 

F3A 

12.66 

16.74 

12.43 

za=5  (~ 

SOOmb) 

GSM 

46.63  678 

42.52 

671 

35.20 

663 

INIT 

38.64 

33.62 

23.71 

ANAL 

36.67 

21.02 

20.69 

F3A 

20.14 

32.73 

19.76 

! 

II 

<0 

N 

200rab) 

GSM 

73.56  639 

68.61 

634 

62.46* 

627 

INIT 

63.39 

58.29 

51.12 

ANAL 

50.32 

35.60 

35.41 

F3A 

34.78 

44.01 

33.88 

Vo=2 

860mb) 

GSM 

6.76(m  s-1)  740 

6.55 

724 

6.50 

709 

INIT 

5.15 

4.82 

4.77 

ANAL 

4.77 

4.39 

4.35 

F3A 

4.05 

4.06 

3.93 

-» 

< 

Q 

II 

X 

575mb) 

GSM 

7.65  705 

7.34 

683 

7.31 

672 

INIT 

6.83 

6.42 

6.39 

ANAL 

6.70 

fa.  25 

6.21 

F3A 

4.39 

4.40 

4.27 

Vo=8  (~ 

225mb) 

GSM 

12.05  627 

11.35 

610 

11.28 

604 

INIT 

10.61 

9.76 

9.70 

ANAL 

10.76 

9.82 

9.77 

F3A 

6.52 

6.54 

6.34 

*rNIT-ANAL  >  GSM-IMIT 


ANALYSIS  TIMES  (GMT) 


Fiq.  3(g).  Same  as  Fig.  3(a)  except  q  at  Layer  1  (~960  mb) 
(top) ,  q  at  Layer  3  (~725  mb)  (middle) ,  and  q  at  Layer  5 
(~435  mb)  (bottom) 


J _ I _ I _ I _ l_i _ L_J 

9/00  10/00  11/00  12/00 

ANALYSIS  TIMES  (GMT) 


ANAL 

INIT 


F3A 


Fig.  3(a)  except  V  at  Layer  B  (~-225  mb) 


two  possible  paths  to  follow.  One  can  either  perform  a  forecast  with  the 
initialized  field,  thereby  ending  the  cycle,  or  one  can  regard  the  ini¬ 
tialized  field  as  a  new  first  guess  field  and  begin  the  cycle  anew.  The 
rationale  behind  the  second  path  is  that  by  repeated  analyses  and  initial¬ 
izations,  one  can  expect  that  if  convergence  occurs,  the  final  fields  will  be 

both  close  to  the  observations  and  in  a  balanced  state.  Indeed,  Williamson 
18 

and  Daley  (1983)  p'"oposed  sequential  iterations  of  01  and  initialization 
for  this  very  reason.  In  our  case,  with  the  model  normal  modes  employed  in 
both  the  analysis  and  initialization,  one  would  hope  that  convergence  would 
be  rapid.  Because  of  computational  constraints,  we  were  restricted  to  no 
more  than  tnree  iterations  per  cycle. 

To  prepare  for  the  iterations,  the  initialized  data  in  the  form  of 
spectral  coefficients  of  divergence,  vorticity,  temperature,  surface  pres¬ 
sure,  and  specific  humidity  had  to  be  converted  into  residuals  of  velocity 
and  the  compound  function  at  the  observation  sites.  This  involved  expanding 
the  spectral  coefficients  to  obtain  velocities  and  heights  at  the  observation 
sites,  subtracting  them  from  the  observed  data,  and  then  repeating  the  analy¬ 
sis  cycle.  When  the  requisite  number  of  iterations  is  completed,  a  forecast 
is  produced  immediately  after  the  initialization  procedure. 

D.  Surface  Pressure  Calculat’ons 

Residuals  are  calculated  at  the  beginning  of  the  cycle  by  means  of  the 
method  outlined  in  Section  II  of  this  report.  That  is,  the  observed  vari¬ 
ables  are  interpolated  by  means  of  01  vertical  correlation  functions  to  the 
a  layers  of  the  model.  The  first  guess  field  is  then  evaluated  at  the 

observation  points  and  subtracted  from  the  observations.  However,  whereas 

surface  pressure  is  analyzed  separately  in  Section  II  of  this  report,  here  it 
is  derived  from  the  heights  near  the  surface  by  a  quadratic  interpolation. 

The  reason  for  this  is  the  excessive  computer  time  needed  to  perform  the  01 
surface  pressure  analysis,  given  the  vast  amount  of  surface  data  available. 
The  interpolation  also  insures  that  the  heights  and  surface  pressure  will  be 


18.  Williamson,  D.  L.  and  R.  Daley,  1983;  A  unified  analysis-initialization 
technique.  Mon .  Wea .  Rev . .  Ill .  1517-1536. 
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in  dynamic  balance,  whereas  an  independent  derivation  of  surface  pressure  may 
lead  to  an  inconsistency  which  would  have  to  be  corrected  by  the  initiali¬ 
zation. 

After  the  analysis  procedure,  the  residuals  of  the  compound  function  P 
are  separated  into  height  and  surface  pressure  residuals.  The  height  resi¬ 
duals  are  then  added  to  the  first  guess  field  to  form  the  updated  heights. 
Surface  pressure  could  have  been  updated  the  same  way  (i.e.,  by  adding  the 
analyzed  residuals  to  the  first  guess  surface  pressure)  but  we  chose  to  in¬ 
terpolate  the  new  heights  to  the  sun  ace  to  form  the  updated  surface  pressure 
field.  This  is  in  keeping  with  NMC's  practice  where  all  surface  pressures 
for  modeling  purposes  are  prescribed  by  the  upper  air  height  observations. 

The  a  surfaces  are  redefined  with  respect  to  pressure  by  the  new 
surface  pressure  before  initialization.  As  in  the  ASAP  program  outlined 
earlier  in  this  report,  the  variables  carried  at  the  a  layers  are  inter¬ 
polated  in  pressure  coordinates  to  the  new  a  layers  with  the  assumption 
that  they  are  linear  in  In  p  (except  for  heights,  which  are  redefined  as 
layer  temperatures  before  interpolation).  With  the  variables  defined  on  the 
new  a  layers,  initialization  can  proceed. 

E.  Run  Stream  Summary 

A  review  of  the  programs  used  in  the  analysis  procedure  may  be  helpful 
in  the  understanding  of  the  full  cycle.  There  are  12  individual  programs  in 
the  current  procedure; 

1.  The  global  spectral  model.  The  model  produces  a  12h  forecast  from 
a  previously  analyzed  field  to  serve  as  the  first  guess  field.  Run  time: 
approximately  290  CPU  seconds  on  the  Cray. 

2.  The  post-processor.  The  spectral  values  of  vorticity,  divergence, 
temperature,  surface  pressure,  and  moisture  are  reproduced  as  velocity,  height, 
surface  pressure,  and  specific  humidity  on  grid  space  on  the  a  layers.  Run 
time;  19  CPU  seconds. 

3.  The  ASAP  residual  calculations.  01  analysis  interpolates  obser¬ 
vations  to  a  layers,  expands  the  spectral  coefficients  of  the  first  guess 
field  at  the  observation  sites,  and  calculates  the  residuals  of  velocities  and 
temperatures.  Run  time;  345  CPU  seconds. 
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4.  Reformation  of  the  residuals.  Missing  observations  are  filled  by 
residuals  equal  to  zero  so  that  all  levels  at  an  observation  site  will  be 
filled.  Also,  grid  boxes  containing  no  observations  are  filled  with  zeroes 
at  their  respective  four  comers  as  in  HT.  Pressure  residuals  are  inter¬ 
polated  from  height  residuals  and  residuals  of  the  compound  function  com¬ 
puted.  Run  time:  63  CPU  seconds. 

5.  Vertical  modes  program.  At  each  observation  site,  the  residuals 
are  projected  onto  eight  of  the  model's  vertical  modes.  Run  time:  4  CPU 
seconds . 

6.  Analysis  program.  The  projected  residuals  are  analyzed  hori¬ 
zontally  by  sequential  fits  to  the  normal  mode  functions.  The  results  are 
the  coefficients  of  the  normal  mode  coefficients  for  all  eight  vertical 
modes.  Run  time:  650  CPU  seconds. 

7.  Evaluation  of  analysis  on  grids.  The  residuals  are  computed  at 
grid  points  by  expansion  of  the  normal  mode  functions  at  the  appropriate 
latitude  and  longitude.  Run  time:  68  CPU  seconds. 

8.  Analyzed  data  on  grid.  Residuals  are  added  to  the  first  guess 
values  computed  at  the  grid  points.  The  compound  function  residuals  are 
divided  into  surface  pressure  and  height  residuals.  First  guess  heights  are 
added  to  each  layer  and  the  updated  heights  interpolated  at  the  surface  to 
produce  surface  pressure  and  the  new  o  structure.  The  heights  are  con¬ 
verted  to  layer  temperature  by  the  Flattery  method.  Run  time:  75  CPU 
seconds . 

9.  The  expansion  of  grid  data.  The  velocities,  temperatures,  and 
surface  pressure  are  converted  to  spectral  coefficients  of  divergence, 
vorticity,  temperature,  and  surface  pressure.  Spectral  coefficients  of 
moisture  are  added  from  the  FGGE  IIl-A  analyses,  while  the  fixed  spectral 
coefficients  of  terrain  height  are  added,  as  well.  Run  time:  11  CPU  seconds. 

10.  Initialization.  The  NMC  initialization  program  is  run,  resulting 
in  spectral  coefficients  of  vorticity,  divergence,  temperature,  surface 
pressure,  and  specific  humidity.  Run  time:  32  CPU  seconds.  If  no  iterations 
are  desired,  or  if  the  full  number  of  iterations  has  been  completed,  the  next 
step  is  a  return  to  the  first  program  where  a  new  12h  forecast  is  produced. 

If  the  data  are  to  be  matched  to  the  initialization  by  iteration,  the  follow¬ 
ing  two  steps  are  added: 
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11.  The  post-processor.  As  with  program  2,  the  initialized  field  is 
converted  to  grid  point  space  in  terms  of  velocities  and  heights  for  use  as 
the  new  first  guess  field.  Run  time:  8  CPU  seconds. 

12.  Residual  recalculation.  The  observations  are  now  matched  with 
the  new  first  guess  field  which  is  expanded  at  the  observation  sites.  Val¬ 
ues  of  zero  are  added  as  before  in  data  void  areas  and  at  missing  layers. 
Surface  pressure  and  height  residuals  are  combined  as  compound  variable  re¬ 
siduals.  Run  time:  380  CPU  seconds.  From  here,  one  returns  to  program  5 
to  repeat  the  cycle. 

F.  Results 

Table  6  compares  analyses  from  the  non-iterative  experiments  with 
model  forecasts  and  with  the  analyses  after  initialization.  It  is  the 
counterpart  to  Table  2  of  Section  II,  except  that  wind  components  are 
separated,  all  12  levels  are  included,  and  only  four  dates  are  available. 

The  large  errors  at  the  highest  levels  led  us  to  believe  that  extrapolation 
by  the  Flattery  scheme,  as  indicated  in  Section  II,  p.  11,  was  a  principal 
source  of  error  in  the  analysis.  Indeed,  the  large  error  seems  to  propagate 
downward  over  the  course  of  the  four  analysis  periods.  In  a  subsequent 
experiment  where  the  upper  level  temperatures  were  replaced  by  first  guess 
temperatures  rather  than  extrapolated  temperatures,  errors  at  the  top  were 
smaller,  but  errors  in  lower  layers  were  unaffected  for  the  most  part. 

There  are  other  important  conclusions  one  may  draw  from  Table  6  when 
contrasted  with  Table  2.  First,  the  differences  between  the  forecasts  and 
the  analyses  are  almost  always  greater  here  than  they  are  in  Table  2.  One 
could  view  this  as  proof  that  our  analysis  has  a  greater  effect  on  the  first 
guess  field  than  an  01  analysis  does.  On  the  other  hand,  it  could  be  an 
indication  of  poorer  forecasts  resulting  from  this  analysis.  Second,  the 
change  produced  by  the  initialization  is  relatively  less  for  this  analysis 
than  for  the  01  analysis,  especially  for  surface  pressure  and  heights  of  the 
upper  o  levels.  This  can  be  understood  in  light  of  the  01  treatment  of 
surface  pressure  and  the  use  of  nmf  for  our  analysis.  In  Section  II  the 
surface  pressure  analysis  was  described  in  detail  as  a  completely  inde¬ 
pendent  analysis,  as  opposed  to  our  derivation  of  pressure  as  an  inter¬ 
polation  (extrapolation)  from  heights  at  the  lower  a  layers.  This  means 
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Table  6 .  RMS  Differences  at  Grid  Points  between  Forecast  and  Analysis  (F-A) ,  Forecast  and 
Initialized  Field  (F-I),  and  Initialized  Field  and  Analysis  (I-A)  for  Heights  (Z) 
at  All  12  Levels  and  for  the  Two  Components  of  Velocity  (U,V)  at  All  12  Layers 
(Surface  Pressure  (PS)  for  the  Four  Analysis  Times  Also  Appears  at  the  Bottom) 
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that  the  gross  atmospheric  mass  structure  as  depicted  by  the  upper  air 
analysis  will  not  necessarily  be  in  dynamic  balance  with  the  surface  analy¬ 
sis  of  Section  II.  In  an  attempt  to  balance  the  atmospheric  structure,  the 
initialization  will  quickly  change  the  surface  pressure  to  be  in  harmony  with 
the  upper  air  mass  structure.  Indeed,  the  initialization  procedure  combines 
heights  and  pressures  into  one  variable  before  adjusting  the  gravity  waves. 
Because  the  overall  structure  of  the  atmosphere  is  mostly  defined  by  the  first 
guess,  with  the  01  performing  some  minor  adjustments,  the  initialization 
procedure  will  substantially  alter  the  pressure  analysis  imposed  by  the 
surface  01  and  substitute  a  field  close  to  the  first  guess.  With  our  analy¬ 
sis,  changes  are  made  to  the  large  scale  structure  by  expansion  in  nmf.  Thus, 
when  NLNMI  takes  place,  there  is  no  tendency  to  draw  the  analysis  back  to  the 
first  guess  field. 

Figs.  4  and  5  are  counterparts  to  Figs.  1  and  2  from  Section  II.  Shown 
are  the  rms  errors  of  forecast,  analysis,  or  initialized  field  against 
observations  for  the  five  dates  of  this  study  (six  dates  if  one  includes  the 
final  forecast  for  1/10/79  OOZ)  at  layers  2,  4,  and  8  for  the  east-west 
component  of  velocity  and  levels  2,  4,  and  8  for  heights.  Vertical  axis 
resolution  varies  from  figure  to  figure  to  ease  comparisons.  Figs.  4  (a-f) 
are  for  observations  retained  after  checking  with  FGGE  III-A  values  as  a  basis 
for  verification.  Figs.  5  (a-f)  use  the  analysis  as  a  basis.  Many  of  the 
figures  indicate  that  the  forecast  sometimes  is  closer  to  the  observations 
than  the  analysis  or  the  initialized  field  for  the  same  time.  This  is 
particularly  true  at  the  higher  levels  and  more  clearly  for  velocities  than 
heights.  At  the  upper  levels,  this  may  be  caused  by  contamination  from  the 
highest  level  due  to  the  Flattery  scheme.  At  lower  levels,  the  error  growth 
is  almost  as  large  as  if  no  update  occurred  at  all.  For  comparison.  Table  4 
of  Section  II  displays  48h  forecast  errors  verifying  at  1/11/79  OOZ.  Taking 
into  consideration  the  differences  in  levels  displayed  for  heights  and  the 
vector  wind  error  versus  the  component  wind,  one  finds  that  the  model  forecast 
errors  are  fairly  close  to  the  error  incurred  by  our  assimilation  cycle.  This 
would  seem  to  indicate  that  not  enough  weight  is  being  given  to  the  observed 
values.  Thus,  despite  the  fact  that  the  analyses  do  somewhat  alter  the  first 
guess  fields  as  is  evident  from  the  figures,  their  impact  is  not  sufficient  to 
draw  the  fields  closer  to  the  observations.  Although  it  was  thought  benefi¬ 
cial  to  filter  the  observations  through  expansion  in  nmf,  it  is  probable 
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Fig.  4(a).  RMS  Differences  as  a  Function  of  Time  between 
the  Analysis  (long  dash).  Initialized  Data  (short  dash)  or 
Forecast  (solid  line)  and  Observations  with  FGGE  III-A  as 
a  Background  Criterion  for  Retention,  for  Z  at  Level  2 
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4(e).  Same  as  Fig.  4(a)  except  for  the  East-West  Compo- 
of  V  at  Layer  4 
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Fig.  7(a).  Difference  Contours  between  the  Forecast  and  Analyzed  Height  Field  for  1/8/79 
OOZ  at  Level  4  (Contour  Interval  is  60  m.  Labels  Are  in  Hectameters) 


Fig.  6(a).  Difference  Contours  between  the  Forecast  and  FGGE  III-A  Height  Field  for  1/8/79 
OOZ  at  Level  4  (Contour  Interval  is  60  m.  Labels  Are  in  Hectameters) 


situation,  and  the  forecast  based  on  the  last  iteration  was  far  inferior  to 
the  forecast  based  on  only  one  iteration. 

There  are  several  possible  explanations  for  the  non-convergence  of 
the  iterative  process.  The  simplest  is  that  the  procedure  is  by  nature 
divergent.  This  would  result  from  either  the  analysis  or  the  initialization 
furthering  the  solution  from  the  data  or  slow  manifold  rather  than  bringing 
it  closer.  In  an  earlier  test  with  just  the  external  vertical  mode  of  Hough 
functions,  this  did  not  seem  to  occur.  In  fact,  the  process  was  seen  to 
converge.  In  our  current  procedure,  however,  the  analysis  and  initializa¬ 
tion  did  not  appreciably  draw  the  first  guess  closer  to  the  observations,  as 
is  obvious  from  Fig.  5,  where,  in  some  cases,  the  analysis  and  initializa¬ 
tion  actually  widen  the  gap  between  the  first  guess  and  the  observations. 
When  this  procedure  was  repeated,  it  merely  reinforced  the  departure  from 
observations.  It  is  also  possible  that  the  introduction  of  higher  fre¬ 
quencies,  just  as  with  the  Machenhauer  scheme  in  the  initialization  pro¬ 
cedure,  can  cause  a  divergence  in  the  solution.  Another  obvious  source  of 
error  is  the  Flattery  scheme  and  its  attendant  extrapolation  at  the  top  of 
the  atmosphere.  But  this  error  was  restricted  to  the  uppermost  levels,  as 
was  discovered  when  the  uppermost  temperatures  were  replaced  by  first  guess 
temperatures  and  divergence  still  occurred.  There  is  still  the  distinct 
possibility  that  a  coding  error  exists,  especially  in  the  recalculation  of 
the  residuals  in  program  12  of  Subsection  E  above.  A  few  extra  tests  are 
necessary  in  order  to  locate  the  source  of  the  growing  errors  and  to  correct 
them.  These  have  not  been  possible  to  date  because  of  budget  considerations 

Figs.  6  (a-e)  show  the  contoured  differences  between  the  five  fore¬ 
cast  heights  and  the  FGGE  III-A  height  fields  at  level  4  for  1/8/79  OOZ, 
l/8/;9  12Z,  1/9/79  OOZ,  1/9/79  12Z,  and  1/10/79  OOZ,  respectively.  The 
contour  interval  is  60  m  and  the  labels  are  in  hectameters.  As  can  be  seen 
there  are  several  areas  where  departures  between  the  two  fields  become  quite 
large  during  the  five  sequences.  Noteworthy  are  the  departures  near  the  US 
East  Coast,  the  US  Rockies,  and  Eastern  Siberia.  The  differences  do  not 
grow  continuously  and  seem  to  drop  off  during  the  last  forecast  period. 

Figs.  7  (a-d)  are  contours  of  differences  between  the  forecast 
heights  at  level  4  and  our  analyzed  heights  at  the  same  level  from  1/8/79 
OOZ  through  1/9/79  12Z,  respectively.  These  may  be  viewed  as  adjustments 
made  to  the  first  guess  field  by  our  analysis,  rather  than  as  a  verification 
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that  too  much  filtering  occurred  and  the  influence  of  the  observations  was 
lost.  This  may  have  occurred  because:  1.  there  are  too  many  data  void 
areas  where  zero  residuals  were  inserted,  2.  there  were  no  weights  given  to 
the  observations  over  the  first  guess  field,  or  3.  the  sequential  calcula¬ 
tion  of  coefficients  has  a  tendency  to  damp  the  wave  amplitudes  resulting  in 
oversmoothing,  as  mentioned  by  HT.  In  fact,  the  lack  of  impact  may  be  a 
result  of  a  combination  of  all  three  reasons  and  further  experimentation  is 
required  to  determine  and  correct  the  precise  cause  of  the  problem.  In¬ 
deed,  the  iterative  experiments  could  have  rectified  some  of  the  problems 
but  they,  too,  failed. 

The  main  goal  of  the  iterative  scheme  is  to  draw  the  analysis  closer 

to  the  observations  while  maintaining  the  large  scale  aspects  of  the 

analysis.  The  iterations  are  akin  to  those  proposed  by  Williamson  and  Daley 
18 

(1983)  ,  who  demonstrated  the  effects  of  repeating  01  analyses  emd  an 

initialization  procedure.  The  expectation  is  that  by  alternating  between  an 
01  which  maintains  geostrophic  balance  and  an  initialization  which  corrects 
the  gravity  wave  components,  an  equilibrium  will  be  reached  such  that  the 
final  analysis  will  be  in  balance  with  respect  to  the  model  and  yet  close  to 
the  observed  data.  Our  attempt  to  perform  this  type  of  Iteration  is  based 
on  the  same  reasoning,  except  that  instead  of  an  01  analysis,  our  analysis 
scheme,  using  only  the  large  scale  nmf  as  interpolation  media,  is  sub¬ 
stituted.  The  large  scale  nmf  would  keep  the  Rossby  modes  and  lower 
frequency  gravity  modes  intact  while  approximating  the  data.  The  initial¬ 
ization  procedure  would  then  adjust  the  higher  frequency  gravity  modes  to 
return  the  fields  to  balance  on  the  "slow”  manifold.  Eventually  the 
iteration  would  hopefully  converge  on  the  ideal  intersection  of  the  "data" 
manifold  and  the  "slow"  manifold  as  depicted  by  Williamson  and  Daley 
(1983).  Unfortunately,  our  iterative  cycle  did  not  converge  at  all.  Table 
7  compares  rms  errors  of  the  analyzed  fields  versus  observations  for  heights 
and  the  east-west  component  of  velocity  for  iterations  2  and  3  with  FGCB 
III-A  data  as  a  basis  for  the  checking  procedure.  The  iterations  were  per¬ 
formed  for  the  1/8/79  OOZ  fields  and  can  be  compared  with  the  first  iter¬ 
ation  depicted  in  Figs.  4  (a-f).  As  is  readily  seen,  there  is  a  definite 
divergence  of  the  iterative  process  with  observed  values  departing  from  the 
analyzed  values.  Subsequent  initialization  did  not  do  much  to  improve  the 
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7/12  8/00  8/12  9/00  9/12 

ANALYSIS  TIMES  (GMT) 

Fig.  5(c).  Same  as  Fig.  4(c)  except  the  Analysis  Is  Used 
as  the  Basis  for  Retention  of  Observations 


ANALYSIS  TIMES  (GMT) 

Fig.  5(b).  Same  as  Fig.  4(b)  except  the  Analysis  Is  Used 
as  the  Basis  for  Retention  of  Observations 
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Fig.  7(b).  Diffi 
12Z  at  Level  4  (' 


Fig.  7(c).  Difference  Contours  between  the  Forecast  and  Analyzed  Height  Field  for  1/9/79 
OOZ  at  Level  4  (Contour  Interval  is  60  m.  Labels  Are  in  Hectameters) 


of  the  forecast.  As  such,  one  would  have  hoped  that  the  largest  differences 
would  occur  in  areas  where  the  forecast  and  FGGE  III-A  differ  most,  as 
depicted  in  Figs.  6  (a-e),  assuming  that  FGGE  III-A  fields  are  sensitive  to 
observations.  Unfortunately,  the  adjustments  depicted  by  Figs.  7  (a-d)  are 
not  major  at  any  time,  even  at  the  final  period  (1/9/79  12Z,  Fig.  6d),  com¬ 
pared  with  the  large  differences  that  appear  between  FGGE  III-A  and  the 
forecast  at  that  time.  In  Eastern  US,  for  instance,  where  the  forecast  has 
heights  well  below  the  FGGE  III-A  analysis  with  a  maximum  difference  of  238 
m,  the  analysis  scheme  only  slightly  increases  the  forecast  heights  by  a 
maximum  of  47.8  m.  Areas  such  as  Eastern  Siberia  or  Antarctica  are  also 
slightly  adjusted  despite  the  large  differences  present  between  the  forecast 
and  FGGE  III-A  data.  Whether  this  pattern  would  continue  indefinitely  is 
unknown,  but  as  residuals  increase,  the  adjustments  to  the  first  guess  field 
should  increase  as  well,  perhaps  halting  the  present  pattern  of  insensi¬ 
tivity  to  observations  and  increasing  forecast  errors. 

G.  Conclusions 

The  current  study  is  not  yet  complete.  Mot  enough  cycles  have  been 
generated  to  test  the  effects  of  the  analysis  and  forecast  cycle  over  the 
long  term.  Nor  has  there  been  enough  testing  to  fully  determine  whether 
even  the  complex  computer  codes  are  genuinely  error  free.  There  are  some 
identifiable  sources  of  error  that  should  be  addressed.  They  are:  1.  error 
due  to  the  analysis  procedure,  2.  vertical  and  horizontal  interpolation 
error,  and  3.  the  Flattery  scheme. 

1.  The  analysis  scheme  employs  an  approximation  that  does  not 
guarantee  a  minimum  error.  As  mentioned  by  HT,  the  error  depends  on  data 
distribution  and  cannot  be  easily  determined  a  priori .  at  present.  Whether 
an  efficient  mathematical  evaluation  of  the  error  can  be  produced  can  only 
be  determined  after  extensive  research  into  linear  algebra  and  statistics. 
The  sequential  procedure  also  tends  to  oversmooth  the  analyzed  fields,  as 
does  the  introduction  of  first  guess  information  (zero  residuals)  in  areas 
that  lack  data.  In  contrast  to  a  least-squares  fit  where  ill-conditioning 
creates  large  uncontrolled  errors  in  data  void  areas,  the  sequential  method 
is  apparently  insensitive  to  data  even  in  areas  of  plentiful  observations. 
Figs.  6  and  7,  in  fact,  indicate  that  the  analysis  scheme  has  very  little 
effect  on  the  first  guess  field,  even  in  the  presence  of  large  observation- 


forecast  differences,  assuming  that  FGGB  III-A  analyses  are  sensitive  to 
observations.  This  seems  to  be  the  casd  even  in  areas  of  considerable  data 
such  as  North  America,  where  one  would  not  expect  zero  residuals  to  be 
retained  at  the  grid  comers  as  is  done  in  the  data  poor  regions.  One 
possible  remedy  for  this,  as  already  mentioned  by  HT,  is  the  introduction  of 
weights  in  the  analysis  procedure.  These  weights  would  be  based  on  statisti¬ 
cal  estimates  of  error  similar  to  01  procedures.  However,  because  nmf  are 
involved  in  estimating  the  observed  data,  the  weights  would  have  to  be 
functions  of  zonal  wavelength  and  frequency.  To  amass  error  statistics  in 
terms  of  wavelength  and  frequency  would  require  reducing  forecast,  analyses, 
and  verifications  into  expansions  of  nmf  and  calculating  the  necessary  means, 
variances,  etc.,  necessary  for  statistical  weights.  This  could  become  a 
formidable  undertaking  unless  some  prior  assumptions  are  made  regarding  the 
nature  of  the  error.  The  weights  could  then  be  modeled,  as  the'^  are  for  01 
techniques,  based  on  the  many  assumptions  and  limited  available  data. 

2.  Errors  due  to  interpolation  are  unavoidable  in  any  analysis 
scheme.  For  our  analysis,  the  chief  sources  of  these  errors  are  the  vertical 
interpolation  with  the  01  scheme  and  the  interpolation  of  the  nmf  to  the 
observation  locations  in  the  horizontal.  If  economy  is  sought,  interpolation 
may  be  substituted  in  the  calculation  of  residuals  at  the  observation  points, 
instead  of  the  current  practice  of  computing  first  guess  fields  directly  from 
the  spectral  expansions.  The  errors  inherent  in  the  vertical  interpolation 
are  discussed  in  the  literature  and  in  the  first  part  of  this  report,  relevant 
to  01  errors.  The  interpolation  of  the  nmf  to  latitudes  of  the  observations 
does  not  necessarily  engender  large  errors,  except  for  the  highest  wave- 
numbers.  It  may  be  possible  to  store  nmf  values  for  the  fixed  locations  of 
upper  air  olservations  and  interpolate  only  to  the  variable  location  sites 
such  as  ship  reports  and  satellite  data.  This  would  necessarily  bias  the 
Northern  Hemisphere  because  of  its  preponderance  of  dry  land  and  developed 
nations  where  most  fixed  observing  stations  are  located.  In  any  analysis 
scheme,  however,  the  Southern  Hemisphere  always  suffers  because  of  the  dearth 
of  verifiable  data  there,  and  there  is  no  way  to  judge  whether  any  further 
lack  of  accuracy  will  be  produced  by  this  shortcut. 
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3.  The  Flattery  scheme,  as  implemented  in  our  analysis,  has  a  fatal 
flaw  which  wreaks  havoc  with  the  uppermost  level.  This  is  due  primarily  to 
the  lapse  rate  in  the  stratosphere  which,  when  extrapolated,  results  in  very 
high  temperature  or  height  values  for  the  top  level.  Corrections  to  the 
Flattery  scheme  were  made  by  substituting  first  guess  temperature  at  the 
highest  level.  Results  showed  a  slight  improvement  of  the  analysis  at  upper 
levels  but  no  significant  change  in  lower  levels. 

Results  from  the  iterative  scheme  have  proved  disappointing  because 
the  procedure  diverged,  producing  analyses  that  were  apparently  farther  from 
the  data  and  requiring  greater  adjustments  during  initialization.  A  recom¬ 
mended  remedy  could  be  iterating  the  analysis  scheme  before  initializing, 
which  would  draw  the  fields  closer  to  observations  but  forfeit  the  benefits 
of  bringing  the  analysis  closer  to  the  "slow"  manifold  by  virtue  of  the 
initialization. 

In  general,  a  full  evaluation  of  our  analysis  procedure  is  impossible 
until  further  testing  and  experimentation  can  be  completed.  To  make  the 
scheme  viable  greater  efficiencies  have  to  be  implemented.  If  vector 
processing  machines  such  as  the  Cray  are  to  be  assigned  the  task  of 
producing  the  analysis,  the  code  must  be  modified  to  allow  full  vector- 
ization.  It  may  also  be  possible  to  economize  by  sticking  to  spectral 
expansions.  The  nmf  are  determined  by  expanding  their  spectral  coeffi¬ 
cients.  It  may  be  possible  to  rewrite  the  analysis  in  terms  of  expansions 
in  Legendre  functions  directly  rather  than  in  nmf  followed  by  re-expansion 
on  a  physical  grid  in  order  to  define  the  proper  spectral  coefficients. 

Also,  it  would  probably  be  beneficial  to  go  directly  from  the  analysis 
procedure  to  the  initialization  without  involving  a  physical  grid.  This 
would  avoid  the  problem  of  determining  surface  pressure  r<  siduals  from  the 
analysis  by  separation  from  the  compound  variable.  This  c<  i  only  be  accom¬ 
plished,  however,  if  the  first  guess  field  is  expanded  in  nmf,  and  the 
velocity  nmf  transformed  to  vorticity  and  divergence  nmf  to  be  compatible 
with  the  initialization.  Despite  the  apparent  magnitude  of  the  challenge,  a 
scheme  that  can  combine  objective  analysis  with  initialization  while  of¬ 
fering  computational  economy  is  a  goal  worthy  of  pursuing. 
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APPENDIX  A.  FLATTERY  ALGORITHM 


Given  the  system  of  linear  constraints  imposed  on  the  interface  temperatures 
{Tj^},  k  =  1,  ....  K  +  1  and  the  layer  temperatures  k  =  1,  ....  K.  We 

assume: 


(1)  that  the  layer  temperature  T.  is  the  arithmetic  average  of  the 


A  A 


bounding  interface  temperatures  Tj^, 


1  A  A 

Tk  =  t  +  Tk+i) 


k  =  1,  . . . ,  K; 


(A-] 


(2)  that  temperature  in  the  layer  bounded  by  two  pressure  levels  p.  , 

A  ^ 


is  linear  in  t.n  p,  so  that  the  interface  temperature  is  given 


in  terms  of  the  layer  temperatures  Tj^  and  as 


A 

T  =  T  W  +  T  W 

k+1  k  L,  k+1  ^k+1  U,  k+1 


where 


W. 


In  (^k+l/^k+l) 


L,k+1 


W. 


U,k+1 


In  (Pk/^k-H) 


(A-: 


k  =  1, 


K; 


(3)  that  both  Eq.  (A-1)  and  oq.  (A-2)  are  valid  also  at  the  lowest 

1  “■* 


A  A 

and  highest  bounding  interfaces,  i.e.,  p.  and  p„. , .  Hence 


A  A 

T  =  2T  -  T 
1  1  2 


2T^  -  (TjWj^  2  +  2^ 


or 


T=(2-W  )T-W  T 

1  L,2^  1  "u,2  2 


(A-3) 


=  T  W  +  T  W 


similarly 


The  constraints  listed  above  may  be  put  into  matrix  form  given  by 


To  solve  Eq.  (A-5)  in  accordance  with  the  least-squares  principle  we  proceed 
as  follows: 

(A)  When  {T^}  is  given, 

A  T  T 

T  =  (A  A)  A  BT  (A-6) 

A 

(B)  When  {Tj^}  is  given, 


T 


T  -1  T  ^ 
(B  B)  B  AT 


(A-7) 


